%% section to load data

%figure out way to read channel name.
clear all


path1='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\Registration and analysis\0. raw single Z\';

path1='C:\Users\Tom�s\Desktop\temp\single Z\';

% get the folder contents
d2 = dir(path1);
% remove all files (isdir property is 0)
dfolders = d2([d2(:).isdir]==1);
% remove '.' and '..'
dfolders = dfolders(~ismember({dfolders(:).name},{'.','..'}));

for s=1:numel(dfolders)
    tic
    f1=waitbar((s-1)/numel(dfolders),['slices completed= ' num2str(s-1) '/' num2str(numel(dfolders))    ])
    path2=dfolders(s).name; % slice name e.g. 's7b3'
    
    dir_to_search=[path1 '\' path2 '\' ];
    
    
    % get the folder contents
    d = dir(dir_to_search);
    % remove '.' and '..'
    d = d(~ismember({d(:).name},{'.','..'}));
    nZplanes=round(numel(d)/3); %there is an extra image for background, so round takes care of that.
%     nZplanes=11
    imageparams=imfinfo([d(1).folder '\'  d(1).name]); % I don't need this since metadata has no info about number of slices and channels for some reason. or at least I can't find it with this function.
    
    %due to memory constraints need to run one channel at a time.
    channels={'Nxph1','Slc6a5','Aldh1a3'};
%         channels={'Sorcs3'};

    slices=path2;
    
    f=waitbar(0,'starting')
    
    
    
    for c=1:numel(channels)
        
        for j=1:nZplanes %loop through z planes
            clear Zmax rawmedfilt raw xind yind masked logicalPeaks mask mask01
            close all
            fclose all
            
            
            fracprogress=((c-1)/3)+ j/(3*nZplanes);
            waitbar(fracprogress, f, {['processing channel ' num2str(c) '/' num2str(numel(channels)) ': ' channels{c} ', Z plane ' num2str(j) '/' num2str(nZplanes)];['Remaining time: ' num2str((1-fracprogress)*(toc/fracprogress)/60) ' minutes']} );
            
            
            %create data matrix for each channel.
            raw=imread([d(j).folder '\' channels{c} '_' num2str(j) '.tif']);
            
            rawmedfilt=raw-medfilt2(raw,[20 20]);
            thresh=10*mad(rawmedfilt(:),0);
            p=FastPeakFind(rawmedfilt,500);
            
            figure;
            imagesc(rawmedfilt); hold on
            plot(p(1:2:end),p(2:2:end),'r.','markersize',8)
            y1=colormap([cmap('magenta',100,0,50)]);
            caxis([200 1000])
%             
%             
          
            
            yind=p(1:2:end); % there is some weird transposition.
            xind=p(2:2:end);
            clear pks pks2
            
            for i=1:size(xind)
                pks(i)=raw(xind(i),yind(i));
                pks2(i)=rawmedfilt(xind(i),yind(i));
            end
            
            SNRraw=mean(pks(:)/mean(raw(:)))
            SNRmedfilt=mean(pks2(:)/mean(rawmedfilt(:)))
            
            
         %% make plot for checking and suplementary figure 2.2   
         %plot raw and medfilt1 side by side.
%             figure
% ax1=subaxis(1,2,1)
% 
% imagesc(raw);
% ax2=subaxis(1,2,2)
% imagesc(rawmedfilt)
% axis square
% y1=colormap([cmap('magenta',100,0,50)]);
% 
% % title('after filtering'
% axis square
% y1=colormap([cmap('magenta',100,0,50)]);
% 
% linkaxes([ax1 ax2], 'xy')


            
            %% create a circular mask radius 3 px
            radius=10; %in pixels
            imageSizeX = 20;
            imageSizeY = 20;
            [columnsInImage, rowsInImage] = meshgrid(1:imageSizeX, 1:imageSizeY);
            % Next create the circle in the image.
            centerX = imageSizeX / 2; % in the middle.
            centerY = imageSizeY / 2;
            circlePixelsLogical = (rowsInImage - centerY).^2 ...
                + (columnsInImage - centerX).^2 <= radius.^2;
            % figure
            % image(circlePixelsLogical) ;
            % colormap([0 0 0; 1 1 1]);
            % title('Binary image of a circle');
            % axis square;
            
            
            logicalPeaks=zeros(size(rawmedfilt));
            for i=1:numel(xind)
                logicalPeaks(xind(i),yind(i))=1;
            end
            
            
            %test convolution, constrain size of operation to the central part so
            %that the size of the mask is the same as the input image. In this case,
            %logicalPeaks is the same size as the original image.
            
            mask=conv2(logicalPeaks,circlePixelsLogical, 'same');
           
            mask01=double(mask>0); %use it for multiplication to apply mask.
            clear mask
            
            %make masked image based on detected puncta
            masked=mask01.*double(rawmedfilt);
%             figure;
%             imagesc(masked);
%             y1=colormap([cmap('magenta',100,0,50)]);
%             caxis([200 1000])
%             
            clear mask01 logicalPeaks rawmedfilt
            
            if j==1
            currentZmax=masked;
            clear masked
            else 
                
            currentmaskedMat=cat(3,masked,currentZmax); % this takes waay too much memory
            currentZmax=max(currentmaskedMat,[],3);
            clear masked currentmaskedMat
            end
            
            
            
            
            
            
            
            
            
            
            
            
            
            % fclose all
            % Slc6a5RawMat(:,:,i)=imread([d(i).folder '\' 'Slc6a5_' num2str(i) '.tif']);
            % fclose all
            % Aldh1a3RawMat(:,:,i)=imread([d(i).folder '\' 'Aldh1a3_' num2str(i) '.tif']);
            % fclose all
        end
        
        %
        Zmax=currentZmax;
        % figure;
        % imagesc(Zmax);
        % y1=colormap([cmap('cyan',100,0,50)]);
        % caxis([200 1000])
        % toc
        clear currentmaskedMat
        
        path3=['C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\Registration and analysis\1.Processed not registered\' path2 '\'];
        
%         imwrite(uint16(Zmax),[path3 channels{c} 'Matlabout.tiff'])
        
        
        
    end
end
waitbar(1, f,'Done!')
waitbar(1, f1,'Done!')
toc









